The EEAaq package allows users to retrieve air quality data for multiple geographical zones, pollutants, and time periods in a single request. Queries are submitted as lists, which enables flexibility in specifying combinations of parameters.
library(EEAaq)
library(tidyverse)
## Warning: il pacchetto 'ggplot2' è stato creato con R versione 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Below we demonstrate the use of query by using different combinations of user-defined arguments.
LAU_IDdata_lau <- EEAaq::EEAaq_get_data(
zone_name = "15146", # LAU zone code
NUTS_level = "LAU", # NUTS level
LAU_ISO = "IT", # Country code for Italy
pollutants = "PM10", # Pollutant
from = "2022-01-01", # Start date
to = "2023-12-31", # End date
verbose = FALSE # Print detailed progress
)
# Preview the first few rows of the dataset
head(data_lau)
LATN_NAME# Identify the names of the areas from which to download the data
zones <- c("Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest","Vlaams Gewest","West-Nederland","Zuid-Nederland")
# Download the corresponding data
data <- EEAaq_get_data(
zone_name = zones, # LAU zone code
NUTS_level = "NUTS1", # NUTS level
pollutants = c("NO2", "PM10"), # Pollutant
from = "2023-01-01", # Start date
to = "2023-12-31", # End date
verbose = FALSE # Print detailed progress
)
unique(data$AirQualityStationEoICode)
## [1] "BELAL01" "BELAT83" "BELHB23" "BETB001" "BETB004" "BETB006" "BETB008"
## [8] "BETB011" "BETBUL1" "BETCHA1" "BETE013" "BETE714" "BETE716" "BETM802"
## [15] "BETMEU1" "BETN043" "BETR001" "BETR002" "BETR012" "BETR701" "BETR702"
## [22] "BETR721" "BETR740" "BETR801" "BETR802" "BETR803" "BETR804" "BETR805"
## [29] "BETR806" "BETR817" "BETR818" "BETR822" "BETR831" "BETR842" "BETR891"
## [36] "BETR897" "BETREG1" "BETVBX1" "BETVBX2" "BETVBX3" "NL00136" "NL00138"
## [43] "NL00236" "NL00237" "NL00240" "NL00241" "NL00247" "NL00546" "NL00551"
## [50] "NL00553" "NL00556" "NL00570" "NL00572" "NL00573" "NL00701" "NL00704"
Note 1: If the query’s zone_name parameter corresponds
to a valid CITY_NAME (i.e., not NULL in the dataset), the
function will return the corresponding data. If no valid CITY_NAME is
associated with the zone_name, the function attempts to retrieve all
available data for the entire country and subsequently filter for the
specified zone_name.
Note 2: For very small towns or certain countries
such as Turkey or Albania, data may not currently be available in the
dataset.This limitation reflects the data unavailability at the EEA
Air Quality Viewer.
Note 3: If the parameters used in the query include
polygon or quadrant, the function outputs an
EEAaq_df_sfc object. Otherwise, it returns an EEAaq_df object, which is
a tibble dataframe.
EEAaq_map_stations generates a static or dynamic map of
user-defined monitoring stations. The function accepts as input either
an object of the EEAaq_df class (default output of the
EEAaq_get_data function), or all other parameters
specifying the area and the pollutants.
EEAaq_df object, the dataset
concerning NO\(_2\) and PM\(_{10}\) in Belgium and The NetherlandsEEAaq_map_stations(
data = data,
bounds_level = "NUTS3",
color = FALSE,
dynamic = FALSE
)
## Simple feature collection with 4 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2.546088 ymin: 50.688 xmax: 6.225231 ymax: 53.18511
## Geodetic CRS: WGS 84
## NUTS_ID LEVL_CODE CNTR_CODE
## 1 BE1 1 BE
## 2 BE2 1 BE
## 3 NL3 1 NL
## 4 NL4 1 NL
## NAME_LATN
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest
## 2 Vlaams Gewest
## 3 West-Nederland
## 4 Zuid-Nederland
## NUTS_NAME MOUNT_TYPE
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest NA
## 2 Vlaams Gewest NA
## 3 West-Nederland NA
## 4 Zuid-Nederland NA
## URBN_TYPE COAST_TYPE geometry
## 1 NA NA MULTIPOLYGON (((4.415738 50...
## 2 NA NA MULTIPOLYGON (((5.776583 50...
## 3 NA NA MULTIPOLYGON (((5.171192 52...
## 4 NA NA MULTIPOLYGON (((5.518671 51...
## points Country ISO AirQualityStationEoICode AirQualityStationNatCode AirQualityStationName Altitude NUTS1 NUTS1_ID NUTS2 NUTS2_ID NUTS3 NUTS3_ID LAU_NAME LAU_ID AirPollutant geometry n
Note: Using the parameter
bounds_level = "NUTS3", the map is generated with internal
boundaries corresponding to the NUTS-3 level. The same output could be
obtained specifying explicitly the zone information.
EEAaq_map_stations(
zone_name = zones,
NUTS_level = "NUTS1",
pollutant = c("NO2", "PM10"),
bounds_level = "NUTS3",
color = FALSE,
dynamic = FALSE
)
## Simple feature collection with 4 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2.546088 ymin: 50.688 xmax: 6.225231 ymax: 53.18511
## Geodetic CRS: WGS 84
## NUTS_ID LEVL_CODE CNTR_CODE
## 1 BE1 1 BE
## 2 BE2 1 BE
## 3 NL3 1 NL
## 4 NL4 1 NL
## NAME_LATN
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest
## 2 Vlaams Gewest
## 3 West-Nederland
## 4 Zuid-Nederland
## NUTS_NAME MOUNT_TYPE
## 1 Région de Bruxelles-Capitale/Brussels Hoofdstedelijk Gewest NA
## 2 Vlaams Gewest NA
## 3 West-Nederland NA
## 4 Zuid-Nederland NA
## URBN_TYPE COAST_TYPE geometry
## 1 NA NA MULTIPOLYGON (((4.415738 50...
## 2 NA NA MULTIPOLYGON (((5.776583 50...
## 3 NA NA MULTIPOLYGON (((5.171192 52...
## 4 NA NA MULTIPOLYGON (((5.518671 51...
## points Country ISO AirQualityStationEoICode AirQualityStationNatCode AirQualityStationName Altitude NUTS1 NUTS1_ID NUTS2 NUTS2_ID NUTS3 NUTS3_ID LAU_NAME LAU_ID AirPollutant geometry n
EEAaq_map_stations(
data = data_lau,
color = TRUE,
dynamic = TRUE,
ID = TRUE
)
## Simple feature collection with 1 feature and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 9.042173 ymin: 45.38829 xmax: 9.277104 ymax: 45.53535
## Geodetic CRS: WGS 84
## GISCO_ID ISO LAU_ID LAU_NAME POP_2021 POP_DENS_2 AREA_KM2 YEAR NUTS3_ID
## 1 IT_015146 IT 15146 Milano 1397715 NA 181.8155 2021 ITC4C
## geometry
## 1 MULTIPOLYGON (((9.120405 45...
## points Country ISO AirQualityStationEoICode AirQualityStationNatCode AirQualityStationName Altitude NUTS1 NUTS1_ID NUTS2 NUTS2_ID NUTS3 NUTS3_ID LAU_NAME LAU_ID AirPollutant geometry n
This function aims to describe the dataset that has been previously
imported, both at a global level, which means considering the complete
set of time stamps and monitoring stations in the dataset, and at the
station-specific level, where summary statistics and information are
grouped by monitoring station.
In addition to basic exploratory descriptive statistics (e.g., average
pollutant concentration, variability, measures of skewness and
kurtosis), the function provides information about the gap length and
the correlation between pollutants if at least two pollutants are
considered simultaneously.
The EEAaq_summary function receives as input an
EEAaq_df object, i.e. the output of the EEAaq get data
function.
summ <- EEAaq_summary(data = data)
## The dataset contains:
## ** 477510 total observations
## ** 56 stations
## ** 8736 time stamps: from 2023-01-01 01:00:00 to 2023-12-31
summ$Summary
summ$Summary_byStat$Mean_byStat
summ$Corr_Matrix
Recall that most pollutants are monitored by EEA on a hourly or daily
basis, posing challenges for interpretation and representation. The
EEAaq_time_aggregate function simplifies this by
aggregating data into annual, monthly, weekly, daily, or hourly
intervals, generating summary statistics for each station in an
EEAaq_taggr_df object.
t_aggr <- EEAaq_time_aggregate(
data = data,
frequency = "monthly",
aggr_fun = c("min", "max", "mean", "median" )
)
t_aggr$TimeAggr
t_aggr$TimeAggr_byPollutant$PM10
To enable quick and intuitive visual analysis, the
EEAaq_idw_map function provides spatial interpolation maps
using the Inverse Distance Weighting (IDW) method (Shepard, 1968). This
technique estimates the value of a variable at unknown locations by
calculating a weighted average of known values, with weights inversely
proportional to the distance from known points. Closer points contribute
more heavily to the estimate, making it a practical approach for
interpolating geolocated air quality data.
EEAaq::EEAaq_idw_map(
data = t_aggr,
pollutant = "PM10",
aggr_fun = "mean",
distinct = TRUE,
gradient = TRUE,
idp = 2
)
## Map initialization started at 2025-02-05 13:30:08.525219
## Map initialization ended at 2025-02-05 13:30:19.841176
## Computing IDW interpolation started at 2025-02-05 13:30:19.841352
## Computing IDW interpolation for: 2023-01-01, 1 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-02-01, 2 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-03-01, 3 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-04-01, 4 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-05-01, 5 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-06-01, 6 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-07-01, 7 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-08-01, 8 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-09-01, 9 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-10-01, 10 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-11-01, 11 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-12-01, 12 of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2025-02-05 13:31:54.439807
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
EEAaq::EEAaq_idw_map(
data = t_aggr$TimeAggr_byPollutant$PM10 %>% dplyr::filter(Date %in% c("2023-01-01","2023-02-01")),
pollutant = "PM10",
aggr_fun = "max",
distinct = TRUE,
gradient = TRUE,
idp = 2
)
## Map initialization started at 2025-02-05 13:32:01.657422
## Map initialization ended at 2025-02-05 13:32:14.605266
## Computing IDW interpolation started at 2025-02-05 13:32:14.605443
## Computing IDW interpolation for: 2023-01-01, 1 of 2
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-02-01, 2 of 2
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2025-02-05 13:32:29.484387
## [[1]]
##
## [[2]]
EEAaq::EEAaq_idw_map(
data = t_aggr,
pollutant = "PM10",
aggr_fun = "mean",
distinct = TRUE,
gradient = FALSE,
fill_NUTS_level = "NUTS3",
dynamic = TRUE
)
## Map initialization started at 2025-02-05 13:32:30.684601
## Map initialization ended at 2025-02-05 13:32:42.167963
## Computing IDW interpolation started at 2025-02-05 13:32:42.16821
## Computing IDW interpolation for: 2023-01-01, 1 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-02-01, 2 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-03-01, 3 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-04-01, 4 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-05-01, 5 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-06-01, 6 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-07-01, 7 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-08-01, 8 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-09-01, 9 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-10-01, 10 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-11-01, 11 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation for: 2023-12-01, 12 out of 12
## [inverse distance weighted interpolation]
## Computing IDW interpolation ended at 2025-02-05 13:34:23.277381
Note: By setting
fill_NUTS_level = "NUTS3", we require the
predictions/interpolations to be aggregated (averaged) at the NUTS-3
level.